clear;
n=1;

rho=0.05;
mu_H=0.4;
mu_L=0.3;
M=1;

up_bound_phi_L=-(mu_L+mu_H)/rho*M
cont=1;
phi_L=-0.1;
k=1;
while cont==1
phi_L=phi_L-0.005;
phi_H=(M-mu_L*phi_L/(mu_L+mu_H))/mu_H*(mu_L+mu_H);

eta_1=(rho+mu_H)/phi_L;
eta_2=(rho+mu_L)/phi_H;

A=[eta_2 -mu_L/phi_H;-mu_H/phi_L eta_1];
om=eig(A);


om_1=om(2);
om_2=om(1);

x_m(k)=max(1/(om_2-om_1)*log(om_1^2/om_2^2*(eta_2-om_2)/(eta_2-om_1)),0);

if x_m(k)<x_m(1)
    cont=0;
end
 phi(k)=phi_L

cons1=(mu_H+mu_L)/(phi_L*mu_L+phi_H*mu_H);
cons2=-(phi_L*mu_L+phi_H*mu_H)/(phi_L*phi_H);
time_to_death_analytical(k) =1/mu_L +phi_H/mu_L*cons1*(exp(cons2*x_m(k))-1);


k=k+1;
end


figure
FSsize=14;
plot([sqrt(mu_L)/sqrt(mu_H)*abs(phi-M)]./M,[ 0 cumsum(time_to_death_analytical(2:end))/sum(time_to_death_analytical(2:end))],[sqrt(mu_L)/sqrt(mu_H)*abs(phi-M)]./M,cumsum(ones(size(phi)))./sum(ones(size(phi))),'r','Linewidth',3)
box('off')
xlabel('coefficient of variation, $\theta$','Interpreter','latex','FontSize',FSsize)
ylabel('cumulative distribution','Interpreter','latex','FontSize',FSsize)
ylim([0 1]) 

annotation('textarrow',[0.5 0.55],[0.6 0.65], 'String','entering firms','Interpreter','latex','FontSize',FSsize)
annotation('textarrow',[0.7 0.75],[0.8 0.85], 'String','stationary distribution','Interpreter','latex','FontSize',FSsize)

